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In 2005, Nagler and Claussen (Phys. Rev. E 71, 067103 (2005)) investigated the time series of the elementary cellular 
automata (ECA) for possible (multi)fractal behavior. They eliminated the polynomial background at b through the direct fit- 
ting of the polynomial coefficients a and b. We here reconsider their work eliminating the polynomial trend by means of the 
multifractal-based detrended fluctuation analysis (MF-DFA) in which the wavelet multiresolution property is employed to filter 
out the trend in a more speedy way than the direct polynomial fitting and also with respect to the wavelet transform modulus 
maxima (WTMM) procedure. In the algorithm, the discrete fast wavelet transform is used to calculate the trend as a local 
feature that enters the so-called details signal. We illustrate our result for three representative ECA rules: 90, 105, and 150. 
We confirm their multifractal behavior and provide our results for the scaling parameters. 

Introduction. - At the present time, a number of different algorithms are well established to analyse 
the singular behavior that may be hidden in time series data, such as the WTMM method [TJ [5J [31 [H [S] , the 
structure function method [5] , the DFA [B] and its variants [JJ [HI H] ■ DFA is a method used to analyze the 
behaviour of the average fluctuations of the data at different scales after removing the local trends. In 2002, 
Kantelhardt et al. [7j provided a generalization of DFA to the case of multifractal time series. Subsequently, 
the latter method started to be widely employed in the literature under the name of MF-DFA. Kantelhardt 
wrote a recent review of the techniques used in processing the fractal and MF time series [5] . 

On the other hand, as already mentioned, a lot of research has been done on fractal signals and ob- 
jects with wavelet transforms (WTs) because the multiscale decompositions implied by the WTs are well 
adapted to evaluate typical self-similarity properties. The efficiency of WTs as 'mathematical microscopes' 
for capturing the local scaling properties of fractals have been noticed since more than two decades [10] . 

It is thus no wonder that there are current efforts towards merging the WTs with DFA procedures [TT] 
as a natural union of powerful tools for quantifying the scaling properties of the fluctuations. In this short 
note, based on this unifying standpoint, which we call WMF-DFA, we focus on the MF properties of ECA 
with periodic boundary conditions. There is only one previous work dedicated to the MF features of ECA 
[12] but there the analysis is performed on the time series of random walk processes generated by some of 
the ECA evolution rules and not directly to the ECA time series as we do here. In addition, Nagler and 
Claussen [13] mention in the final part of their work the possibility of considering their spectral analysis for 
MF signals instead of monofractal ones. We recall that many important applications of ECA are in biology, 
chemistry, and soft materials, where MF properties are to be expected. For example, an interpretation of 
ECA rules 90 and 150 can be made in the context of catalytic processes [13], also the rule 126 can be used 
as a conceptual model of biological cell growth [T3]. On the other hand, the rule 110 is interesting because 
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it has been proven that any mathematical algorithm can be mapped to a ECA having this rule. It is also 
considered as an intrinsic generator of randomness. 

WMF-DFA. - An important advantage of the vanishing moment property of wavelets (see the Appendix) 
is that it helps detrending the data. We are interested in revealing the MF properties [TS] of ECA. To 
separate the trend from fluctuations in the ECA time series, we follow the discrete wavelet method proposed 
by Manimaran et al. [TT]. This method exploits the fact that the low-pass version resembles the original 
data in an "averaged" manner in different resolutions. Instead of a polynomial fit, we consider the different 
versions of the low-pass coefficients to calculate the "local" trend. Let x(tk) be a time series type of data, 
where = kAt and k — 1,2, ... , N. Then the algorithm that we employ contains the following steps (for 
more mathematical details, see the Appendix): 

1. Determine the profile Y(k) — X)i=i( 2; (^i) — ( x )) °f the time series, which is the cumulative sum of the 
series from which the series mean value is subtracted. 

2. Compute the fast wavelet transform (FWT), i.e., the multilevel wavelet decomposition of the profile. 
For each level m, we get the fluctuations of the Y(k) by subtracting the "local" trend of the Y data, 
i.e., AY(k; m) = Y(k) — Y(k; m), where Y(k; m) is the reconstructed profile after removal of successive 
details coefficients at each level m. These fluctuations at level m are subdivided into windows, i.e., into 
M s = int(N/s) non-overlapping segments of length s. This division is performed starting from both 
the beginning and the end of the fluctuations series (i.e., one has 2M S segments). Next, one calculates 
the local variances associated to each window v 

F 2 (v,s;m) = varAF((^ - l)s + j;m) , j = l,...,s, v = l,...,2M s , M s =mt(N/s). (1) 

3. Calculate a qth order fluctuation function defined as 
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Fq{3 ' m) = { 2M S S \ F2 ^ s; m )l ?/2 ] ( 2 ) 

where q G Z with g / 0. Because of the diverging exponent when q — > we employed in this limit a 

f 1 2Ms 1 
logarithmic averaging F$(s; m) — exp I — — Vln|F 2 (i/,s;m) } as in [JJ [TH] . 

I 2Ms v=i J 

In order to determine if the analysed time series have a fractal scaling behavior, the fluctuation function 
F q (s;m) should reveal a power law scaling 

F q (s;m) ~ s h ^\ (3) 

where h(q) is called the generalized Hurst exponent [16j since it can depend on q, while the original Hurst 
exponent is h(2). If h is constant for all q then the time series is monofractal, otherwise it has a MF behavior. 
In the latter case, one can calculate various other MF scaling exponents, such as r(q) and /(a) [15] . 

Application to ECA. - We apply the previous algorithm to the time series of three illustrative ECA 
as classified by Wolfram in 1984 17J. The chosen rules are the following: 90, 105, and 150. For the first and 
the last rules the updates are given by 



n+1 ]mod2, (4) 

where r = and r = 1, respectively. It is well known that rule 90 has the appearance of a Sierpinski triangle 
when responding to an impulse (first row is all 0s with a 1 in the center) . Nagler and Claussen [13] found that 
the rule 150 displays a Sierpinski-like self-similar structure of fractal dimension (If = 1.69 (golden mean) 
instead of the standard one of 1.58. In a subsequent paper [18] . Claussen showed that its time behavior can 
be solved as a two-step vectorial, or string, iteration, which can be viewed as a generalization of Fibonacci 
iteration generating the time series from a sequence of vectors of increasing length. This could explain the 
difference in the fractal dimension. As for the rule 105, it is known to be complementary to the rule 150, 
i.e., /io5 = 1 — /i5o, where / is the neighborhood-depending updating rule. 
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Fig. 1: Rule 90: (a) Time series of the row signal. Only the first 2 7 points are shown of the whole set of 2 17 data 
points, (b) Prohle Y of the row signal, (c) Generalized Hurst exponent h(q). (d) The r exponent, r(q) — qh(q) — 1. 
(e) The singularity spectrum /(a) = q — T (<l)- The calculations of the multifractal quantities h, r, and /(a) are 
performed both with MF-DFA and the wavelet-based WMF-DFA. 



Table 1: The Hurst exponent h(2) for the three ECA rules examined in this paper as obtained by means of the 
MF-DFA and WMF-DFA methods and in each of the cases for four initial (first row) pulses as indicated. 





MF-DFA 


WMF-DFA 




I 

(■ ■ -010 ■ ■ ■ ) 


II 

(■ ■ ■ 0110 ■ ■ ■ ) 


III 

(■ ■ ■ 01010 ■ ■ ■ ) 


IV 

(■ ■ ■ OHIO ■ ■ ■ ) 


I 

(• ■ ■ 010 ■ ■ ■ ) 


II 

(■ ■ ■ 0110 ■ ■ ■ ) 


III 

(■ ■ ■ 01010 ■ ■ ■ ) 


IV 

(■ ■ ■ OHIO ■ ■ ■ ) 


R90 
R150 
R105 


0.8972 
0.9427 
0.9427 


0.8972 
0.9413 
0.9413 


0.8898 
0.9541 
0.9542 


0.9451 
0.9296 
0.9296 


0.8961 
0.9293 
0.9294 


0.8961 
0.9529 
0.9529 


0.9229 
0.9513 
0.9514 


0.9787 
0.9407 
0.9407 



We have analysed the time series of the so-called row sum (or total activity) ECA signals, i.e., the sum 
of ones in sequences of rows, employing Daubechies (Db) wavelets. Various types of Db wavelets have been 
used but we have found that a better matching of the results given by the WMF-DFA method with those 
of other methods is provided by the Db-4 wavelets with four filter coefficients. Our results are illustrated 
in Figs. (|TJ)— d3j) - The fact that the generalized Hurst exponent is not a constant horizontal line is indicative 
of a multifractal behavior in all three cases. In addition, the fact that the r index is not of a single slope 
is another clear feature of multifractality. The values of the Hurst exponent h(2) for four types of initial 
conditions are given in Table 1. We also present the corresponding fluctuation function F2 in Fig. Q for 
the impulsive initial condition. The strength of the multifractality is roughly measured with the width 
Aa = cv max — cv m j n of the parabolic singularity spectrum /(a) on the a axis. For example, for the impulsive 
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Fig. 2: Same plots as in Fig. [I] but for rule 150. 



initial condition, Aa 90 = 0.99 98(1.0 1 32), Aai 50 = 1.011(1.00 75), and Aai 05 = 1.00 83(1.03 2 5) when the 
MF-DFA (WMF-DFA) are employed. We notice that the most "frequent" singularity for all the analysed 
time series occurs at a = 0.568, where the width Aa of rule 90 is shifted to the right with respect to those 
of 105 and 150. According to our results, the strongest singularity, a m i n , of all time series corresponds to 
the rule 90 and the weakest singularity, a ma x, to the rule 150. 

In conclusion, in general terms, our algorithm implementation shows that embedding the discrete wavelet 
transform in the MF-DFA technique is a well-suited procedure to analyze the multifractal properties of the 
ECA. Indeed, we get similar results to the other methods but computationally faster because we employ a 
lesser number of windows. Our results represent a confirmation of the fact that ECA patterns of different 
magnitudes follow different scaling laws, i.e., the ECA have intrinsic multifractality that does not depend 
on the set of initial data that we used. Therefore, when processes thought to be multifractal are simulated 
with (E)CA, their intrinsic multifractal behavior should be taken into account as a feature of the simulation 
procedure rather than of the multifractal behavior of the simulated processes. 
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Fig. 3: Same plots as in the previous figures but for rule 105. 



Appendix 

WT: continuous and discrete. - The WT of a function or distribution function x(t) is given by 

W x {a,b) = - x{t)^[ )dt, (5) 

where if is the analyzing wavelet, b € M. is a translation parameter, whereas a 6 M + (a ^ 0) is a dilation 
or scale parameter, and the bar symbol denotes complex conjugation. One fundamental property that we 
require in order to analyze the singular behavior of a signal is that ip(t) has enough vanishing moments [51ll9j. 
A wavelet has n vanishing moments if and only if it satisfies /_ t k ip(t)dt = for k = 0, 1, . . . ,n — 1 and 
I-oo t k ip{t)dt ^ for k = n. This means that a wavelet with n vanishing moments is orthogonal to all 
polynomials up to order n — 1. Thus, the WT of x(t) performed with a wavelet ip(t) with n vanishing 
moments is nothing else but a "smoothed version" of the n'th derivative of x(t) on various scales. 

Since the ECA data are notoriously discrete it is important to consider a discrete version of ([5]) . Generally 
the orthogonal (discrete) wavelet transform (DWT) is employed. This is only one of the different forms of 
WTs 201, by which the wavelets are associated to orthonormal bases of L 2 (M). In this case, the wavelet 
transform is performed only on a discrete grid of the parameters of dilation and translation, i.e., a and b take 
on only integer values. In fact, for the numerical implementation of the DWT the multiresolution analysis 
(MRA) has been introduced. 

The representation of a function or process x{t) with the DWT is given in terms of shifted and dilated 
versions of the wavelet function ip(t), and its associated scaling function <p(t) (20j [2T] . Within this framework 
and considering that the scaling and wavelet functions 

l p m>n (t) = 2 m / 2 <p(2 m t-n) ) ^ m , n (t)=2 m ^(2 m t~n), m,neZ (6) 
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Fig. 4: Log-log plot of the fluctuation function F2 versus scale for: (a) rule 90, (b) rule 150, and (c) rule 105. 



form an orthonormal basis, then one can write the expansion of x(t) as follows 



M-l 



= X! a m ,n¥rn ,n{t) + ^ d m , n l/j m . n (t) , (7) 



n \ m-rno / 

where the scaling or approximation coefficients a mn , and the wavelet coefficients d m ^ n are defined as 



x(t)ip m . n (t)dt, d m , n = J x(t)^ m ^ n (t)dt, (8) 
with m and n denoting the dilation and translation indices, respectively. 



FWT. - To calculate a TO) „ and d m ,m Mallat [5D] developed the FWT in which the MRA approach is 
involved. The FWT algorithm connects, in an elegant way, wavelets and filter banks, where the multires- 
olution signal decomposition of a signal X, based on successive decomposition, is composed by a series of 
approximations and details which become increasingly coarse. At the beginning, the signal is split into an 
approximation and a detail part that together yield the original. The subdivision is such that the approxi- 
mation signal contains the low frequencies, while the detail signal collects the remaining high frequencies. By 
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repeated application of this subdivision rule on the approximation, details of increasingly coarse resolution 
are separated out, while the approximation itself grows coarser and coarser. 

The FWT calculates the scaling and wavelet coefficients at scale m from the scaling coefficients at the 
next finer scale m+1 using the following formulas 

a m , n = /] h[k - 2n]a m+ i,fc, (9) 

k 

d m ,n = g[k - 2n]a m +i,k, (10) 

k 

where h[n] and g[n] are typically called low pass and high pass filters in the associated analysis filter bank. 
In fact, the signals et m ,„ and d m ^ n are the convolutions of a m+ i „ with the filters h[n] and g[n] followed by a 
downsampling of factor 2 [20], respectively. 

Conversely, a reconstruction of the original scaling coefficients a m +i jn can be made from the following 
combination of the scaling and wavelet coefficients at a coarse scale 

Qw+i,n = / , (h[2k - n]a mik + g[2k - n]d mik ) ■ (11) 

k 

It corresponds to the synthesis filter bank. This part can be viewed as the discrete convolutions between 
the upsampled signal a m j and the filters h[n] and g[n], that is, following an upsampling of factor 2 the 
convolutions between the upsampled signal and the filters h[n] and g[n] are calculated. The number of levels 
depends on the length of the signal, i.e., a signal with 2 L values can be decomposed into (L + 1) levels. To 
initialize the FWT, we consider a discrete time signal X = {x[l], x[2], . . . , x[N]} of length N = 2 L . The first 
application of ([9|) and (fT0|) . beginning with a m+ i „ = x[n], defines the first level of the FWT of X. The 
process goes on, always adopting the m + 1th scaling coefficients to calculate the "m" scaling and wavelet 
coefficients. Iterating ([9]) and (fT0|) M times, the transformed signal consists of M sets of wavelet coefficients 
at scales m = 1, . . . , M, and a signal set of scaling coefficients at scale M. There are exactly 2^ L ~ m ^ wavelet 
coefficients d TO , n at each scale m, and 2^ L ~ M ' scaling coefficients a,M,n- The maximum number of iterations 
is M max = L. A three- level decomposition process of the FWT is shown in Figure [5] 
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Fig. 5: The structure of a three-level FWT. 
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